###########################################################
##
## Distance to transcontinental and important roads
## 
## Use transcon2008 Jedwab AND Storeygard 2020
## distance to "transcon2008" = 1 (transcontinental)
## distance to "transcon2008" = 2 (important)
##
###########################################################

output_version <- '2024-07-22'

## fishnet
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"

fishnet_adm0_only.prj <- fishnet_adm0_only %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.fishnet.zones <- fishnet_adm0_only.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
fishnet.zones.pts = st_as_sf(as.data.frame(M.fishnet.zones), coords = c("X", "Y"), crs = 4326)
fishnet.zones = st_drop_geometry(subset(fishnet_adm0_only,select=c(OBJECTID, GID_0)))
fishnet.zones.pts = dplyr::bind_cols(fishnet.zones.pts,fishnet.zones)

##
## Transcontinental 2008
## distance to "transcon2008" = 1 (transcontinental)
## distance to "transcon2008" = 2 (important)
lake_chad_sf_fn <- file.path(wkdir,"local/gis_data/018_transportation/roads_NW_transcon_2008.shp")
if(!file.exists(lake_chad_sf_fn)){
  mapinfo_fn <- file.path(wkdir,"local/gis_data/018_transportation/roads_NW_transcon_2008.TAB")
  mapinfo_data <- st_read(dsn = dirname(mapinfo_fn), layer = "roads_NW_transcon_2008")
  
  mapinfo_sf <- mapinfo_data %>% st_as_sf()
  #mapinfo_lines <- st_collection_extract(mapinfo_sf, "LINESTRING")
  
  table(mapinfo_sf$country)
  table(st_geometry_type(mapinfo_data))
  st_crs(mapinfo_data,4326)
  
  
  lake_chad_sf <- mapinfo_sf %>%
    dplyr::filter(country %in% c('Cameroon','Chad','Niger','Nigeria'))
  
  table(lake_chad_sf$country)
  table(st_geometry_type(lake_chad_sf))
  
  st_write(lake_chad_sf,lake_chad_sf_fn)
} else {
  lake_chad_sf <- st_read(lake_chad_sf_fn)
}

# subset lines as transcontinental
transcon2008 <- lake_chad_sf %>%
  dplyr::filter(trn2008 == 1)
st_crs(transcon2008) <- st_crs(fishnet.zones.pts)

fishnet.zones.pts$dist2transcon2008<-NA
fishnet.zones.pts$dist2transcon2008 <- unlist(nngeo::st_nn(fishnet.zones.pts, transcon2008, k = 1, returnDist = T)$dist) / 1000

# subset lines as important
important2008 <- lake_chad_sf %>%
  dplyr::filter(trn2008 == 2)
st_crs(important2008) <- st_crs(fishnet.zones.pts)

fishnet.zones.pts$dist2important2008<-NA
fishnet.zones.pts$dist2important2008 <- unlist(nngeo::st_nn(fishnet.zones.pts, important2008, k = 1, returnDist = T)$dist) / 1000

## distance to roads in Lake Chad ##
dist2transimp_roads_fn = file.path(wkdir,"proc_data",sprintf("DIST2TRANS_IMP_roads_%s.dta",output_version))
dist2transimp_roads_out <- fishnet.zones.pts %>%
  sf::st_drop_geometry() %>%
  dplyr::select(OBJECTID,GID_0,dist2transcon2008,dist2important2008) %>%
  dplyr::rename_with(.,tolower) 

# Adding the label value
attr(dist2transimp_roads_out$dist2transcon2008, "label") <- "Euclidean distance (km) to transcontinental roads (2008) : Michelin map"
attr(dist2transimp_roads_out$dist2important2008, "label") <- "Euclidean distance (km) to important roads (2008) : Michelin map"
attr(dist2transimp_roads_out, "label") <- sprintf("dist2transcontinental_roads.R last run on %s",Sys.Date())
head(dist2transimp_roads_out)

## Save the data
haven::write_dta(dist2transimp_roads_out, dist2transimp_roads_fn)